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EFFICIENT LOW-SPEED FLIGHT IN A WIND FIELD 


Michael A. Feldman 
(ABSTRACT) 

A new software tool was needed for flight planning of a high altitude, low speed un- 
manned aerial vehicle which would be flying in winds close to the actual airspeed of the 
vehicle. An energy modeled NLP formulation was used to obtain results for a variety of 
missions and wind profiles. The energy constraint derived included terms due to the wind 
field and the performance index was a weighted combination of the amount of fuel used 
and the final time. With no emphasis on time and with no winds the vehicle was found 
to fly at maximum lift to drag velocity, V m d . When flying in tail winds the velocity was 
less than V m d , while flying in head winds the velocity was higher than V m d . A family of 
solutions was found with varying times of flight and varying fuel amounts consumed which 
will aid the operator in choosing a flight plan depending on a desired landing time. At 
certain parts of the flight, the turning terms in the energy constraint equation were found 
to be significant. An analysis of a simpler vertical plane cruise optimal control problem 
was used to explain some of the characteristics of the vertical plane NLP results. 
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Chapter 1 
Introduction 


Problems associated with finding the ‘best’ path for an object have long been of interest to 
both mathematicians and practitioners. The classical brachistochrone problem, proposed 
by John Bernoulli in 1696 [1] is an example of an early problem that led to the development 
of the Calculus of Variations. In more recent times, the problem of optimal ascent of 
a rocket-powered vehicle proposed by Robert Goddard (see [2] for example), provided 
motivation for modern development in trajectory optimization. 

As noted optimal shaping of aerospace trajectories has provided the motivation for much 
recent study of modern optimization theory and algorithms. Current industrial practice 
favors approaches where the continuous-time optimal control problem is transcribed to a 
finite-dimensional Non-Linear Programming problem (NLP) by a discretization process. 
Two such general formulations are implemented in the POST [3] and the OTIS [4] codes. 

There are several existing codes for various flight planning applications. The OPTIM 
[5] code, developed by Analytical Mechanics Associates, Inc., and the MITRE [6] code were 
written to be used for flight planning purposes or as part of an on-board flight manage- 
ment system for general aviation or turbojet transport aircraft. Both codes minimize fuel 
consumption or direct operating costs. Only the vertical plane aspects of the flight are 
optimized in these codes. The horizontal position coordinates that the vehicle follows are 
predetermined depending on the locations of navigational aids and airports. 


1.1 Flight Problem 

The motivation for this work is flight-planning for the Theseus vehicle, built by Aurora 
Flight Sciences. The aircraft is to be used to sample the atmosphere over the South Pole 
and is to operate for the most part in an autonomous fashion. The aircraft data specific 
to Theseus are included in Appendix C. While the numerical results are specific to the 
Theseus airplane, the approach is rather general and applies to similar vehicles which fly 
at low speeds in wind fields. The nominal mission includes take-off, long-range flight to 
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a remote site of interest near the South Pole, a scientific data-collection leg and return 
to base. A ground-based human pilot will perform take off and landing but the aircraft 
acts like a robot during the remainder of its flight. It is controlled by providing specific 
coordinates or waypoints to fly through and by providing speed and altitude commands. 
This waypoint idea will be discussed further in Chapter 3. 

The science mission requires flight over a prescribed path at a prescribed altitude. We 
are not focusing on trajectory shaping for this part of the flight since it is mostly defined by 
the scientific data requirements. The vehicle is to operate at high altitudes (60,000 - 80,000 
feet) during this science leg of the mission. The vehicle will fly from Christchurch, New 
Zealand to the vicinity of the South Pole, take data, and fly back. It is about 2800 nautical 
miles from the base to the Pole. The wind-field is presumed known (and persistent); its 
magnitude is expected to be a significant fraction of the vehicle’s true airspeed. Average 
winds at these altitudes are substantial as will be seen by a figure in Chapter 5. A wind 
model obtained from NASA Goddard Space Flight Center describes horizontal winds which 
are in a clockwise direction around the Pole. The wind speed depends on the distance from 
the pole and on the altitude. 

The specific problem for flight planning is getting back to the base. Given the winds 
and the current ‘ state’ of the vehicle can it return safely? Since the vehicle is landed by a 
human operator it is preferable to land during daylight hours. In addition, since the vehicle 
is flimsy due to its low wing loading, we do not want to land in bad weather. Thus, we 
are interested in the range of possible times when we can get back. Will there be sufficient 
light? Will the local weather be okay? To answer these questions we pose a family of flight 
optimization problems where the performance index is a weighted combination of fuel and 
time (i.e. J = Wf + ji tf,fi > 0). The weighting parameter fi defines a one-parameter 
family of flight optimization problems. There should be a family of solutions of varying 
times to return home and fuel consumed. The limits will correspond to the longest time to 
return home which will also be the most fuel efficient path and the shortest time to return 
home. 

In the present paper a discretization is used that is specially adapted to the flight 
problem of interest. Among the unique aspects of the present discretization are: a least- 
squares formulation for certain kinematic constraints; the use of energy ideas to enforce 
Newton’s Laws; and, the inclusion of large magnitude horizontal winds. In order to offer 
more flexibility in the presence of the strong horizontal winds, the horizontal path is not 
fixed as it was in the OPTIM and MITRE codes. Energy exchange due to turning will also 
be included. 
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1.2 Overview 


In order to solve the flight planning problem a mathematical model is necessary. Various 
levels of fidelity are available such as point-mass, energy, and cruise modeling. In the 
next chapter a description of the various models available and reasons to use each one are 
provided. Following this is the development of the finite dimensional energy model. This 
development includes the derivations of the constraint equations to use in the computer 
code to solve the problem. Chapters 4 and 5 discuss some results from running vertical 
plane and full three dimensional versions of the energy model code. Chapter 6 explains 
some of the results using a simple, vertical-plane, continuous optimal control problem. 
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Chapter 2 
Modeling 


An appropriate mathematical model is needed in solving the current flight planning prob- 
lem. In stability and control a highly detailed mathematical model including rigid body 
mechanics is usually used. However, this is a flight performance problem so only motion of 
the vehicle as a whole is of interest. Positions and motions of the control surfaces on the 
vehicle, for example, are not important in this type of problem so the most detailed model 
needed in flight performance is a point mass model. Simplified, reduced order versions of 
the point mass model are also available and these include the energy model and the cruise 
model. 


2.1 Point Mass Model 

The derivation of the point mass model for symmetric flight over a flat, non-rotating Earth 
with horizontal winds is in Appendix A. While the point-mass model is high fidelity, it 
is not best for our purposes. The reason is that in the point-mass model the controls 
are directly related to forces - small disturbance motions include the phugoid mode and 
four zero eigenvalues (3 corresponding to position and one to heading). Thus the system 
has 4 poles at the origin and is very difficult to control. Many waypoints close together 
would be needed to capture the dynamics of the phugoid oscillations for example. A more 
approximate model with less dynamics is necessary in this case. Here we focus on the 
energy model first developed by Kaiser [7] . 

2.2 Energy Model 

In this model velocity and altitude are combined into a single variable energy. The absolute 
energy is defined to be the sum of the potential and kinetic energy per unit weight. The term 
absolute is to remind us that the kinetic energy is computed from the inertial velocity, Vj, 
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as defined in Appendix A. The differential equations for altitude and path angle evolution 
are equilibrated, reducing the dynamic order of the model. In the vertical plane case five 
state variables (two position, two velocity and weight) are reduced to three state variables 
(range, energy and weight). An energy constraint equation ensures that the change in total 
energy from one point to another is equal to the work done by the net external force acting 
on the vehicle. This model is simpler than the point-mass model and is the one which will 
be used to develop the computer code. 

2.3 Cruise Model 

To get a cruise model the energy model is simplified even more by equilibrating the energy 
state. This leaves only two state variables (range and fuel weight consumed). This model 
will be used for an analysis of the continuous optimal control problem to characterize the 
nature of optimal paths in Chapter 6. 
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Chapter 3 

Finite Dimensional Energy Model 


In keeping with modern practice [3, 4] we shall approximate the flight problem by a finite- 
dimensional transcription. Rather than considering the state/control variables as ‘arbi- 
trary’ functions of time, we shall prescribe their values at certain points, along with ap- 
propriate rules for interpolating between these points (in our case a mixture of piecewise 
linear and piecewise constant approximation). This approach has the effect of replacing 
the infinite-dimensional problem (find the state/control functions) with a finite-dimensional 
problem (find certain state/control values). Since the models are generally nonlinear we 
have a Non-Linear Programming problem (NLP). 

Our goal is to develop an efficient NLP formulation for the flight planning problem, 
including the important effects of winds. The NLP framework requires that we identify 
a parameterization of the flight path, a scalar-valued cost functional and vector-valued 
constraints. The constraint functions must enforce the important physics of the flight 
problem to a reasonable level of approximation. 

It is convenient to parameterize the set of possible flight paths in terms of quantities 
that are useful in the flight guidance problem. As noted previously, the guidance of the 
vehicle is naturally implemented in terms of waypoints. Certain waypoints through which 
the airplane will fly are specified. Figure 3.1 illustrates a horizontal projection of four such 
waypoints. At each waypoint we consider the following variables: 

• x, y, h horizontal position coordinates and altitude 

• V, X speed and heading, relative to the air mass 

• E total mechanical energy per unit weight 

• m remaining mass 

• £ throttle setting 
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Figure 3.1: Projection of Waypoints in the Horizontal Plane 
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For conceptual purposes we consider a point-mass model wherein the dynamics are 
described by Newton’s Laws. The States’ and ‘controls’ will be recorded at each waypoint. 
The variables include position coordinates, velocity with respect to the surrounding air 
mass, heading angle and vehicle mass. All of these variables, except heading angle, are 
imagined to change smoothly between waypoints. The heading angle will be assumed to be 
nearly constant between waypoints; any change will occur in a short period of time in the 
vicinity of each waypoint. For concreteness it is imagined that this heading change occurs 
as we approach a waypoint. For the bulk of the flight from waypoint i to point 2 + 1 the 
relative velocity heading angle is constant at the value x%- Controls include the throttle 
setting, £, and are taken to be piecewise constant from a ‘departing’ waypoint to the next 
(‘arriving’) waypoint and then to change instantly at the ‘arriving’ waypoint. The wind 
velocity is a smooth function of position. 

We envision that the waypoints are well separated so that the flight path inclination 
(7) is small. Vertical equilibrium of the aircraft will require equality of lift and weight. 
Additionally, along this path the net work done on the vehicle must equal the change in 
kinetic energy. Finally, the weight change must balance the integral of the fuel rate. A 
detailed description of the various constraints is now presented. 


3.1 Horizontal Path Constraint 

Here the requirement that the integrated path segment must connect adjacent waypoints is 
enforced. For notational convenience, let r represent the horizontal position components; 
and, let v represent the horizontal plane projection of the vehicle’s velocity (relative to the 
air mass). V w is the wind-field which has been assumed horizontal. 

We have 


r = v + V w or f i+ 1 - n - 


r l *+ 1 


v(t) + V w (r(t ), h(t)) dt = 0 


(3.1) 


The integral in (3.1) is approximated by averaging the integrand over the flight segment. 
For the moment we defer a detailed description of this averaging and simply introduce the 
symbol vf for its value. 

The constraint (3.1) is written as: 


ei = f i+ 1 - fi - vf(t i+ 1 -ti) = 0 


The vector quantities may be computed from the data at the waypoints. Note, however, 
that we have elected not to include time in our parameter list. Instead, we solve for 
A ti = (i i+ 1 — ti) from the least-squares problem 

min ||e*d| 2 . 

At 


8 



This leads to the requirement 


With this A U in hand, one can compute the norm of the residual error from 

||e*|| 2 = ||An|| 2 - ||uf At;|| 2 

where Af^ = rj+i — T{. It is tempting to impose the constraint ||ei|| 2 = 0. Unfortunately, 
since this constraint is quadratic, we find that when the constraint vanishes, so does its 
gradient. This will introduce a loss-of-rank feature in the Kuhn-Tucker optimality system 
at the solution point and cause considerable trouble in the NLP problem (see [8], p 78 
ff). Instead we opt for a different pathology and impose the scalar constraint that the y 
component of the residual e* vanish. This leads to 

(Vi +1 - Vi) \\vf\\ 2 - vfy (Afi, vf) = 0 (3.3) 

The choice of y component over x component is arbitrary. Note that for some geometries 
the constraint (3.3) may be satisfied while the path is not kinematically feasible. We must 
test any numerical ‘solution’ of our problem to ensure this situation does not occur. We 
will also include an inequality constraint 

A U - A t min > 0 (3.4) 

where At m i n is some prescribed minimum allowable time. The constraint (3.4) ensures that 
time goes forward and, loosely, that the waypoints are not ‘too close’. 

3.1.1 Calculation of vf 

We use an implicit Euler rule and calculate the integrand as the simple average of ‘left’ 
and ‘right’ values. That is, we have 

«? = ( 1 / 2 ) (<? + <?), 

with 

vf = V w (fi) + Vi (cos Xi e-x + sinXi e y ) 
vf = V w (f i+l ) + V i+1 (cos Xi e x + sin Xi e y ) ■ 

Note that we have used the heading Xi i n both the vf' and the iff 1 expressions, whereas 
the remaining data are local to the waypoint. Our model supposes that the changes in 
local airspeed and wind-field are smooth, whereas the heading change is concentrated as 
we ‘arrive’ at the waypoint. For all but the last few seconds of this flight segment we have 

x(t) = Xi- 


(3.2) 
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3.2 Energy Change Constraint 

We start with Newton’s Second Law. 

mV/ = F (3.5) 

where Vi is the inertial velocity: 

F/ = V + V w (3.6) 

Distance traveled is expressed by equation (3.6) multiplied by the change in time. The 
force in equation (3.5) is dotted with this distance to get an expression for energy. 

mVj . Vf dt = F • (V + V w ) dt (3.7) 

To continue with the derivation of the energy constraint equation two special cases will be 
considered separately: wings level flight and turning flight. Each of these conditions will 
affect the force in equation (3.7) differently. 

3.2.1 Wings-Level Flight 

We now analyze the wings-level flight segment between adjacent waypoints. Note this is 
the bulk of the flight since the turning is supposed to happen in the brief time preceding 
arrival at the new waypoint. Expanding the force in equation (3.5) yields the following 
equation. 

F = (T - D) e t + L e n + mg e z (3.8) 

where 

e t is a unit vector tangential to F, 

e n is a unit vector normal to V in the vertical plane, 

e z is a unit vector in the direction of gravity. 

These unit vectors are illustrated in Figure 3.2. The drag, D, in equation (3.8) is the level 
flight drag with a load factor of n = 1. Substituting equation (3.8) into equation (3.7): 

mVi • V/ dt= [(T - D)e t + Le n ] • Vdt+[(T - D)e t + Le n ] • V w dt + mge z • (V + V w ) dt (3.9) 
Since V = Ve t , equation (3.9) is simplified slightly. 

mVi • Vi dt = (T — D) V dt+ [(T — D)e t + Le n ] • V w dt + mge z • (V + V w ) dt (3.10) 

The absolute specific energy is defined to be the sum of the potential and kinetic energy 
per unit weight. 

£ = _L(vy^) +/ * (3.ii) 
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Figure 3.2: Unit Vectors 
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Taking the derivative with respect to time, the following is obtained: 


dE 1 .a T -* , 

- = - ( v,.n )+h 

After substituting h = — {V + K«) • 4 and multiplying both sides by the differential dt we 
get 


dE = ■ VA dt — (V + V w ) ■ e z dt 

9 

Combining equations (3.10) and (3.12) with simple manipulations leads to: 

'T-D' 


(3.12) 


dE = 


mg 


V dt 


T-D„ L „ 

H e n 

mg mg 


• V w dt 


This form is integrated to yield 

(T-D) 


J dE =! 


mg 


V dt 


T-D„ L „ 

e t H 

mg mg 


■ V w dt 


(3.13) 


The integrand of the first term on the right side of equation (3.13) is the specific excess 
power, P s = (T — D)V/mg. The product T V is the power available , P a , while the product 
D V is the power required , P r . The other two terms in equation (3.13) are given temporary 
names so they can be evaluated separately. 


Tl = ( t-d 




To = 


mg 

— t -LA e . v 

mgjVn V w , 


W 


(3.14) 


Equation (3.13) becomes 

j dE= [P s dt + j [Ti +T 2 ] dt 

From Figure 3.2 we see that the tangent vector can be represented as 


e t = cos 7 e h - sin7 e z 

where is a unit vector in the direction of the horizontal projection of the aircraft velocity. 
Ti now becomes 

Ti = [ cos 7 (4 • V w ) ~ sin 7 (e z ■ K»)] (3.15) 

However, in our model the winds are in the horizontal direction only, so when the wind 
vector is dotted with e z , a vertical unit vector, the second term in the brackets disappears. 
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The horizontal unit vector, eh, can be written as eh = cosx e x + sin% e y and substituted 
into equation (3.15). 

(T — D\ 

T\ = cos 7 cos (xw - x) V w 

V mg ) 

where Xw is the heading angle of the wind velocity and V w = || || . 7\ is multiplied and 

divided by the aircraft velocity, V , to get an expression in terms of specific excess power 
which can be calculated. 

T\ = (-y^) [ cos 7 cos (Xw - X)\ 

Now the T 2 term will be evaluated. Since the vehicle is assumed to be in symmetric flight, 
lift can be written as 

L = mg cos 7 (3.16) 

The normal unit vector, e n , can be written as 

e n = - sin 7 e h - cos 7 e z (3.17) 

Equations (3.16) and (3.17) are substituted into expression (3.14) to yield: 

T 2 = cos 7 [sin 7(^/1 • Vw) + cos7(e z • V w )]. (3.18) 

Again, because the winds are horizontal, equation (3.18) simplifies to 

T 2 = — cos 7 sin 7 cos(%^ - x)V w - 

This equation is multiplied and divided by the aircraft velocity, V , so that V sin 7 can be 
written as the vertical velocity, h, and T 2 becomes 

T 2 = -~cos^fCos(xw - x)h 

Using the assumption that the path angle (7) is small and cos 7 « 1, T\ and T 2 can be 
written as 

T\ = P s y cos( Xw -x) (3-19) 

r 2 = -^cos (Xw~x)h (3.20) 

We de fin e a new variable, T = ^ cos(x«, — x). After substituting T into equations (3.19) 
and (3.20), T\ and T 2 finally become 

T\ = 

T 2 = 


TPs 

-Th 
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Tih 


(3.21) 


The absolute energy integral can now be written as 


P s dt + P s T dt — 


i iE =f 

Rearranging terms, the energy balance equation is 

J dE = J P s ( 1 + P)dt- J T dh 


(3.22) 


Note that for the case with no winds T = 0 and equation (3.22) reduces to the usual energy 
integral result. 


3.2.2 Turning Segment 

There is also energy exchange in the turn in the vicinity of each waypoint. During the turn 
there will be a horizontal component of the lift vector denoted by \\Lh\\ = Z/sin/q where 
li is the vehicle bank angle. This force is dotted with the inertial velocity and integrated 
over time to get the work done in the turn. 

Work due to banked turn = f Lh • (V + V w )dt 


The force is perpendicular to the vehicle velocity, Lh • V = 0, while the contribution from 
the wind velocity leads to 


Work due to banked turn 


L sin fji V w cos(x + 


7T 

2 


Xw) dt. 


(3.23) 


Instead of integrating over time it is more convenient to integrate over the change in heading 
angle. To produce this change of independent variable the third equation in (A. 23) used. 
Assuming that the winds are constant near the waypoint and that the flight path angle, 7, 
is zero, we observe the familiar result: 


mV x = Tsin ji 


or in differential form: 

L sin fi dt = mV dx (3.24) 

Equation (3.24) is now substituted into equation (3.23) to get an equation which can be 
integrated over the heading angle. 


Work due to banked turn = / mVV w [— sin(x — Xw)} dx 

= mVV w [cos(x w - X)\\ b a- 


(3.25) 


One special case of this work in a banked turn would be similar to a problem presented 
in Halfman [9]. Suppose the airplane starts out flying directly against the wind at the 


14 





Top View 


Figure 3.3: Front View of Airplane in Banked Turn 
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same velocity magnitude as the wind. After a 180° turn the plane would be flying in the 
same direction as the wind and at the same relative velocity magnitude. This means that 
x(b) = Xw, x( a ) — Xw ~ tt , V = V W . Using equation (3.25) in this example, the work due 
to banked turn = 2 mV 2 . This is consistent with the difference in kinetic energy before 
and after the turn. Note that in this example there is no change in potential energy. The 
change in kinetic energy = \m{ 2V) 2 — -m( 0) 2 — 2 mV 2 . 

Up until now the assumption was made that the airplane was in vertical equilibrium and 
therefore the load factor was equal to one. However, there is additional work being done 
during the turn near a waypoint due to a load factor greater than one. This increased load 
factor increases drag so it affects the specific excess power in the vicinity of each waypoint 
where turning occurs. We can create a similar integral to the first term on the right side 
of equation (3.22). This will be an integral of the difference of the power required for the 
increased drag and the power required for the level-flight drag over the time it takes to 
make the turn. 

r At d( 1) — Din) 

Work due to maneuvering drag = / — V[1 + T\ dt (3.26) 

Jo W 

where At is the time to complete the turn and W is the weight. To find this difference in 
drag we start with a representation of drag as 

D = C D (ipF 2 ) 5 

A parabolic drag polar is being used so 

Cd = Cd 0 + kCi (3.27) 


The lift coefficient can be expressed in terms of the load factor (n) as 


C L = 


pV 2 S/2 


From this the drag equation becomes 


D = -pV 2 S C Do + k 


pV 2 S/2 


(3.28) 


Then the drag in equation (3.28) is subtracted from the drag for a load factor of one and 
simplified. 

m~D(n) = ipV*Sk [(;v%) 2 - 

= fS(i-" 2 ) 
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Now the integral from equation (3.26) can be evaluated. 


Work due to maneuvering drag 


r At 2 kW 

L ~pVS 


(1 — n 2 )[ 1 + T\ dt 


We assume the turn is made in a circular path with respect to the surrounding air so the 
load factor will be constant. Everything else in the integral will also be constant during 
the turn except for the T term so that is all that remains in the integral. 


2 kW, 2x 

Work due to maneuvering drag = -(1 — n ) 


After integrating we have the following 

2 kW 


Work due to maneuvering drag = 


pVS 


(l-n 2 ) 


pVS 


At ■ 


At ■ 


rAt 


T dt 


V w (- sin(x w - Xr) + sin(xw - Xi)) 


V 


X 


(3.29) 

where the subscripts, r and Z, refer to the right and left waypoints respectively. While the 
change in time is not known directly, it can be calculated by dividing the change in heading 
angle by the heading rate (%). Assuming a horizontal coordinated turn we have 


X = —Vn 2 — 1 

A y V 


(3.30) 


Substituting equation (3.30) into equation (3.29) 


Work due to maneuvering drag 


lfl(- sin(Xu» - Xr) + sin(x™ - Xi))l] 


(3.31) 


The final energy constraint equation comes from combining the total energy from equa- 
tions (3.22), (3.25) and (3.31) and using an implicit Euler scheme to approximate the 
integrals. 

0 = E r — Ei 


— Ai(l + Ei)(Ps r + Psi)/2 
+ i(h r — hi) 

+ m r Vr V Wr (cos(xr - Xwr) - cos (xi - Xw r )) 


- MI X =T (IXr - »l 

+ ^1 sin(x«v - Xl) - sin(x„, - Xr)|) 


(3.32) 


3.3 Mass Change Constraint 

The final equality constraint is an integral form of equation (A. 24). This requires that the 
change in weight between adjacent waypoints equal to the integral of the fuel-rate. The 
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fuel- rate is dependent on the specifies of the engine. In this case those details are outlined 
in Appendix C. We have 

m r — mi + J Wf dt = 0. 

As before, we use an implicit Euler scheme to approximate the integral in terms of data 
values at the adjacent waypoints. 

m r - mi + ^(W fr +W fl ) = 0. (3.33) 

The flight time is given by (3.2). 
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Chapter 4 

Vertical Plane Flight Numerical 
Results 


A code was written which uses simplified vertical plane versions of the constraints from the 
previous chapter. Several problems were run on this version of the code for a few reasons. 
First of all this is a simpler case and was used to see if the results made sense. Secondly 
head winds and tail winds were introduced to see how they affect the flight path. 


4.1 Description of Problem 


The states are range, cr, energy, E , and mass, m. The controls in the problem are altitude, 
h, and throttle setting, £. Since the energy can be computed from equation (3.11) we 
include the vehicle’s velocity with respect to the surrounding air mass in the analysis. 

The NLP problem was set up in the following manner. The range coordinates were kept 
in a separate array and are not available to the optimization algorithm. It is convenient 
to arrange the remaining flight variables in a matrix. Each column of the matrix contains 
the relevant variables at a waypoint. The variables marked with an * are the optimization 
variables while the rest are fixed by the boundary conditions. 


X = 

[®0 *1 • 

• • x Np ] 

■ h. 

*h 2 

h-N P 

Vi 

*V2 

... v Np 

m i 

*m 2 

• • • *m Np 

. *£i 

*6 

• • • *£n p 


The optimization was subject to energy and mass constraints between every waypoint. 
A simplified version of equation (3.32) is used for the energy constraint. Since the flight 
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is restricted to a vertical plane, the turning section of the energy constraint is not needed. 
Now T = where both V w and V are positive in the positive x direction. 


0 = E r — Ei 

- dd(l + Ti){Ps r + Psi)/2 (4.2) 

+ J~i(h r — hi) 


The mass constraint is the same as equation (3.33). With these two constraints the following 
summarizes the numbers of waypoints, variables, constraints, and degrees of freedom in the 
NLP problem. 


• N p waypoints 

• 4 • (N p — 2) + 3 variables 

• 2 • (N p — 1) equality constraints 

• 2 • N p — 3 degrees of freedom 


For our numerical studies there are 29 waypoints, so there are 111 optimization variables, 
56 constraints, and 55 degrees of freedom. 

The time from one waypoint to the next is calculated by dividing the distance traveled 
by the average inertial velocity at which the aircraft flies. 

Xi+i -Xi- J V h dt = 0. (4.3) 


The integral is approximated as 


/ 


V h dt 


2 


(4.4) 


which is substituted into equation (4.3) and solve for At. The throttle setting, £, must 
remain between 0.1 and 1 while the altitude, h, had a lower bound of 0. 

All the results were run for the case of flying from x = —4000 nmi to x = 0 nmi over 
a fixed grid of points. The initial condition on altitude was 10,000 feet while the initial 
condition on velocity was 150 ft/sec. There are final conditions imposed which bring the 
vehicle back to the initial velocity and altitude. As mentioned in Chapter 1 the performance 
index was a linear combination of the weight of fuel consumed and the total flight time, 
J = Wf + /i tf where fi> 0. 
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4.2 Results 


Figures 4.1 through 4.4 are the results of running the vertical plane version of the code 
with no winds and varying the parameter, /i, in the cost function (J = Wf + n tf). As 
emphasis on time increases, the total time to fly decreases and the plane uses more fuel 
because it has a higher throttle setting for a greater percentage of the flight. In Figure 4.2 
(intermediate fi) and Figure 4.3 (large ji) the throttle stays at 100% even after the airplane 
starts to descend, using both gravity and maximum power for more speed. 

At the fifteenth waypoint on Figure 4.1, h = 76, 400 ft and the weight is 4942 lb. From 
equation (B.3) the minimum drag velocity is found to be 416.9 ft/sec at that altitude and 
weight. From the optimal solution, the velocity at the fifteenth waypoint is actually 416.9 
ft/sec so the airplane is flying at the minimum drag velocity. The data at the rest of the 
waypoints for this case with no winds and no emphasis on time reveal similar results. 

Figure 4.4 shows the family of solutions mentioned in the Introduction. It was con- 
structed by taking final times and final fuel consumed from the cases run while varying 
fi. It should be noted that only the points marked with a ‘plus’ sign in the figure are 
from actual data. The lines drawn between them are not from data but are only there 
to better illustrate the general behavior exhibited by varying the emphasis on time. This 
figure shows that as ji increases from 0 to 2.26 lb/sec, the flight time decreases from about 
22 hours to about 16 hours. The final weight also decreases from 4600 lbs to about 4420 
lbs. The final weight decreases because more fuel is used when the throttle is at 100% for 
a greater portion of the flight as fi increases. 

Figures 4.5 through 4.8 are the results of running the code with a tail wind of 50 ft/sec 
and varying the emphasis on time in the cost function. By applying the same analysis as 
the case with no winds to find the minimum drag velocity, the eighth waypoint on Figure 
4.5 is found to have u = V/V m d = 0.926 so the plane is flying at less than the minimum 
drag velocity. This is characteristic of all of the waypoints in this case with a constant tail 
wind and no emphasis on time. 

Figures 4.9 through 4.12 are the results of running the code with a head wind of 50 
ft / sec and varying the emphasis on time in the cost function. Waypoint number fifteen on 
Figure 4.9 is at u = 1.016 so the vehicle is flying slightly above the velocity for minimum 
drag and this is typical of the entire flight. By looking at Figures 4.4, 4.8, and 4.12, it is 
obvious that as the wind against the airplane increases, the range of possible flight times 
decreases. 

These vertical plane results took from 9 seconds to 55 seconds on a DEC Alpha 2100 
computer system. In general the run time was shorter as the emphasis on time increased. 
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Figure 4.9: Vertical Plane, 
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Figure 4.10: Vertical Plane, Head Wind = 50 ft/sec, fi — 0.0226 lb/sec 















Chapter 5 

Three Dimensional Flight Numerical 
Results 


5.1 Description of Problem 

The states in this three-dimensional problem are horizontal position coordinates, x and 
y, energy, E , and mass, m. The controls in the problem are altitude, h, throttle setting, 

£, and heading, As before, the path is parameterized in terms of flight speed and the 
energy is computed from equation (3.11). 

The path is approximated by specifying at each waypoint the seven values: (x, y, h, V , m, x, £). 
Between each adjacent pair of waypoints we enforce the equality constraints from kinematics 
(3.3); from energy (3.32) and from mass change (3.33). In addition, we have the inequality 
constraint (3.4). The bounds on the controls are the same as those for the vertical plane 
problem in the previous chapter. For a simple flight problem we suppose that the values 
of the first five variables are specified at the initial point and that the final values of the 
horizontal position, velocity and heading angle are specified. The rest of the variables, each 
denoted by an asterick in equation (5.1), are the optimization variables. 

In this three dimensional problem the horizontal coordinates cannot be held constant 
over the entire flight as they were in the vertical plane problem because turning at the 
waypoints must be allowed to occur so all of the variables are conveniently viewed in one 
matrix, A. 


Xl 

*x 2 

• • • Xn p 

y i 

*y 2 

■■■ Vn p 

hi 

*h 2 

hN p 

Vi 

*V 2 

... v Np 

rrii 

*m 2 

■ ■ ■ *m Np 

*Xi 

*X2 

• • • Xn p 

*£i 

*6 

■ ■ ■ *€n p 


34 



with 


• N p waypoints 

• 7 • (N p — 2) + 4 variables 

• 3 • (N p — 1) equality constraints 

• 4 • N p — 7 degrees of freedom 

For our numerical studies there are 15 waypoints so there are 95 optimization variables, 
42 constraints, and 53 degrees of freedom. For this three dimensional problem winds are 
always horizontal in a counterclockwise direction about the origin. The problem will start 
at some distance from the origin at a high altitude with some approximated amount of fuel 
left to simulate the aircraft at the end of the data-taking leg of the mission. The results 
from the previous chapter can be used to approximate how much fuel will be left after 
taking data. In all cases in this chapter 1000 lbs of fuel will be the amount said to have 
been consumed before heading home. All of the constraints are used as they are derived in 
Chapter 3. There is also an inequality constraint (3.4) to keep the waypoints from getting 
too close together. 

The problem was solved for two main cases of initial position as seen in Figure 5.1. 
The first case starts at point A where x = 0 nmi and y = —1000 nmi and the second 
starts at point B with x — 0 nmi and y = 1000 nmi. The final destination is always at 
x = —2800 nmi and y — 0 nmi. The airplane is started at an altitude of 60,000 ft, since it 
will be at a high altitude when completing the science mission. There are final conditions 
imposed which bring the vehicle to an altitude of 10,000 ft and an airspeed of 150 ft/sec. 

The problem is run with two different models. The first is an analytic model with 
wind speeds varying as in Figure 5.2. The winds increase linearly with distance from the 
Pole and piecewise linearly with altitude and they are circulating west to east about the 
Pole. The maximum wind speed is 100 ft/sec which occurs 2800 nmi from Pole above 
40,000 feet. The first few results are with the start of flight at point A and varying /q the 
emphasis on time in the cost function. A second case will start at point B also varying 
fi. The second wind model is a cubic spline fit of wind data obtained from the NASA 
Goddard Space Flight Center (GSFC). In this model the winds are still horizontal only 
and they are also circulating west to east about the Pole. The data is taken from average 
daily values, and wind speed is described as a function of distance from the Pole (latitude) 
and altitude. The wind field from —40° to —90° latitude is pictured in Figure 5.3. Since 
data is only available between altitudes of about 18,300 feet and 102,000 feet, the initial 
altitude in these problems is 60,000 feet as before but the final altitude is 19,000 feet. The 
performance index was the same linear combination of the weight of fuel consumed and the 
total flight time, J = Wf + fi tf as in Chapter 4. 
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Figure 5.1: Horizontal Projection of Three Dimensional Problem 
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Figure 5.2: Winds Varying with Range and Altitude 
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5.2 Results 


Figure 5.4 shows the horizontal projections of the results with fi — 0 and fi = 3.048 lb/sec 
with the flight beginning at point A in Figure 5.1. With this starting point the direction 
of the wind is helping the vehicle return to the base. Other intermediate cases were run 
and are discussed later but only the extreme cases are shown in this figure for clarity. The 
short line vectors coming from each waypoint, denoted by a circle, shows the heading of the 
airplane with respect to the surrounding air. With fi — 0 or a small number the airplane 
starts its flight by almost turning away from the destination to fly a more curved path. 
This maneuver is made to allow the circulating winds to help later in the flight. As /i is 
increased the airplane takes a more direct path to its destination. 

Figures 5.5 through 5.7 show range histories of altitude and time histories of velocity, 
throttle setting, and weight obtained by varying the emphasis on time with the flight 
beginning at point A. With fi = 0 in Figure 5.5 the altitude does not drop below 40,000 ft 
for most of the flight since this is the minimum altitude where the winds remain at their 
maximum for the given range. In the same figure some oscillation in the throttle setting 
is apparent. This oscillation will be discussed in the Summary and Conclusions Chapter. 
These figures show that as ji increases in Figures 5.6 and 5.7, the throttle is near 100% for 
a greater portion of the flight just as for the vertical plane results. 

Figure 5.8 shows the time - fuel performance trade-off. This type of plot can be very 
useful when deciding which path to use depending on time constraints. As noted in the 
Introduction, operational issues are expected to dictate arrival times. If the vehicle needs 
to be back as soon as possible, due to an approaching storm for example, one of the paths 
at the lower left of the plot should be chosen. It should be noted that this figure shows 
that there is a definite limit to how fast the plane can return home. The code can produce 
solutions which use more fuel by increasing the emphasis on time in the performance index 
but the time to return home will not decrease anymore past a certain point. This is because 
the throttle is already at 100% for the entire flight. There will also be conditions where it 
may be possible for the vehicle to get home faster but if it tries it may run out of fuel before 
it gets there. Of course if the operator is not worried about the time of arrival, he should 
choose a path at the right hand side of the plot to minimize fuel used, and so provide a 
fuel margin. 

Figure 5.9 shows the horizontal projections of the results with \i = 0 and /i = 3.048 
lb/sec with the flight beginning at point B in Figure 5.1. Now the circulating winds are 
mainly against the direction of the flight. The top curve in this figure is the one for fi = 0. 
The sharp turn at the third to the last waypoint is where the plane dropped below 10,000 
ft (there are no winds). This is why the heading is in the same direction as the actual 
flight path for the last few waypoints. Figures 5.10 through 5.12 show results with varying 
the emphasis on time with the flight beginning at point B. These figures show that as fi 
increases, the throttle is near 100% for a greater portion of the flight just as for the vertical 
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plane results. Figure 5.13 shows the performance trade-off. This time there is a smaller 
range of times (between 12 and 23 hours) to return the plane to the base. 

These three dimensional results took from 4 seconds to 30 seconds on a DEC Alpha 
2100 computer system. As in the vertical plane cases, in general the run time was shorter 
as the emphasis on time increased. 

Figures 5.14 through 5.18 show results using the NASA GSFC wind model (Figure 5.3) 
starting at point B. In this case the code was run with the wind magnitudes at 50% the 
value in Figure 5.3 because in several of the cases the airplane was not able to return to 
the base with the given amount of fuel in those high winds. 

With this much variability in the wind speeds with altitude and latitude there is the 
possibility of local minima. However, this possibility was not considered in this analysis. 
In the future a variety of initial guesses on the variables can be used to try to explore the 
possibility of local minima. The results are similar to those using the analytic model. 

Figure 5.14 shows that there isn’t much variation in the horizontal projections of the 
paths when fi is varied from 0 to a large number. As in using the analytical wind model, 
Figures 5.15 through 5.17 show the throttle setting staying at its maximum of unity for a 
greater portion of the flight, time decreasing, and final weight decreasing as fi increases. 
Figure 5.18 shows that there is a smaller range of flight times, between 11.75 and 14.5 
hours, in the family of solutions. 
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Chapter 6 

Continuous Optimal Control 
Formulation in the Vertical Plane 


In this chapter a version of the flight optimization problem is studied in the context of 
optimal control. The intent is to find the structure of extremal paths from the optimality 
conditions with the objective of better understanding the numerical results detailed in 
Chapters 4 and 5. The results will not be used to calculate extremal paths; our numerical 
results are based on the NLP formulations described earlier. 


6.1 Optimal Control Formulation of the Cruise Model 

To get to the vertical plane cruise model we begin with the cruise model in Appendix A. 
First it is reduced in order by realizing that there is only one horizontal position coordinate 
now for range, x. Heading angle is always constant so x 1S set equal to zero and the x 
equation is eliminated. The bank angle fi will also be set equal to zero since the aircraft 
will not be turning. Two states are equilibrated: h = 0 implies that there is a small path 
angle or 7 = 0 and 7 means L = W. As mentioned in Section 2.2, velocity and altitude 
are combined into an energy state. 

x = V + V w 

Wf = ^ ( 6 . 1 ) 

a \w w ) 

To get a cruise model the energy model is simplified even more by setting E 
results in P a = P, ■ Now only two states remain: 

x = V + V w 

W f = ^ 

J r\ 

The controls are velocity, V, and altitude, h. 


= 0. This 
( 6 . 2 ) 
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6.1.1 Derivation of the Breguet Range Equation 


The Breguet range equation will now be derived for the case of no winds. In the cruise 
model the independent variable may be changed from time to range (x) so that there is 
only one state equation, 


W>, = Wj_ = cP. 
* dx rjV 


(6.3) 


After integrating equation (6.3) and using a similar analysis to Anderson’s [10] by assuming 
steady, level, unaccelerated flight and that 77, L/D, and c are constant throughout the flight, 
the Breguet range equation is obtained. 


R 


7] L w 0 
— — in — 
c D w\ 


(6.4) 


where wo is the initial weight of the aircraft and w\ is the final weight. This shows that 
the airplane should instantaneously fly at (j$) in order to maximize the range. The 
corresponding velocity for minimum drag is as follows. 


Vmd ~ 


2 W 


\ pS(C 


L) rnd 


(6.5) 


By looking at equation(6.5) it can be seen that as the weight of the aircraft decreases the 
minimum drag velocity will decrease. The velocity is also dependent on the air density so 
as the altitude increases, the density decreases and the velocity increases. Figure 6.1 shows 
these relationships graphically with curves of constant weight on a plot of altitude versus 
minimum drag velocity. 

This analysis begins to explain why the vehicle tends to fly at the minimum drag velocity 
with fi = 0 in the absence of winds as discussed earlier. Also the same case shows in Figure 
4.1 that as the altitude increases the velocity increases which agrees with the previous 
analysis of equation(6.5). 

When head winds and tail winds are added to the problem, basic flight mechanics [11] 
says that the airplane should fly at u < 1 in a tail wind and u > 1 in a head wind. This is 
consistent with what is seen in Section 4.2. 


6.1.2 Hodograph 

The term hodograph was coined by Hamilton in 1845 [12]. He was interested in the way the 
relative velocity vector (represented as a point in a plane) changed as two bodies moved 
under their mutual gravitational attraction. Hamilton observed that under Newtonian 
inverse-square gravitation the graph of the velocity vector is always a circle. In modern 
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control the term hodograph is used to describe the set of possible state-rates as the con- 
trol is allowed to vary over admissible values. Using the ordinary differential equation 
representation, 

Z = 

we are interested in the set of points produced by the function / as u varies at fixed z. For 
the current application the state consists of two components namely x and Wf. 

Getting a picture of the hodograph for the cruise model is a fairly simple task. There 
are two controls, velocity and altitude, so a two parameter family of points can be drawn. 
These are plotted as a family of curves; Each curve is plotted for a constant altitude while 
varying the vehicle airspeed, V, from a lower bound at the stall speed to an upper bound 
which causes the power required to exceed the maximum power available. Each curve 
is similar to the power required versus velocity curve since Wf is proportional to power 
required. The hodograph is pictured in Figure 6.2. The range rate and fuel rate values in 
this figure are nondimensional. The range rate is scaled by the velocity for minimum drag 
at standard sea level and standard weight, V m d- The fuel rate is scaled by the standard 
weight divided by the time scale, W s /T sc \ , where T sc \ = V m d/g- It is drawn for the case 
where there are no winds. The curve which starts at the far left of the figure is the one 
corresponding to the lower bound due to the terrain, in this case h — 0 or sea level. As the 
altitude increases the curves move to the right and get shorter due to the velocity limits 
already mentioned. As shown in the figure there is an envelope created by these curves 
which appears to be linear in the middle portion between the range rates of approximately 
1 and 2. 

In order to show that the envelope is linear, Taylor and Mann’s [13] method for deter- 
mining envelopes of plane curves is used. We have a family of curves in the (</>, ip) plane. 
Each curve is parameterized by a value of a. That is F(<f>, ip, a) = 0 at fixed a determines 
a curve in the plane. In the present case we have x = <p, Wf = ip and h = a. 

x = v 

Wf = ^Pr(V, h) = fi(V, h) (6,6) 

where P r (V,h) is the power required. A derivation of the power required equation is 
included in Appendix B. Rearranging equation (6.6) to be in the form F(<p,ip,a) = 0 
yields the following: 

F(V, Wf , h) = h{V, h)-Wf = 0 (6.7) 

where V is used for x. The condition on the envelope from Taylor and Mann [13] is 

F s (<P,iP,a) = 0 (6.8) 

where F 3 = When this condition is applied to equation (6.7), we find that ^ = 0. 
Solving this equation results in V = V m d • When this value of V is substituted back into 
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equation (6.6), the ratio of ^ (the slope of the envelope) turns out to be a constant for a 
constant weight. This shows that the envelope is a straight line. The curved portions are 
only due to the limits on altitude. So, because V = Vmd , the linear portion of the envelope 
corresponds to the airplane flying at maximum lift to drag. 

6.1.3 Optimal Control Problem 

There are lower and upper bounds on both controls. The lower bound on h will depend on 
the terrain being overflown. The lower bound on V is the stall speed of the airplane which 
is dependent on the altitude. There is an upper bound on the lift coefficient which has an 
effect on the maximum altitude at which the airplane can fly. The specific excess power 
must be greater than or equal to zero. This means that the maximum available power 
must be greater than or equal to the power required, P arnax > P r - The throttle setting, £, 
is multipled by the maximum power available and can be adjusted between 0 and 1 so that 

^^max = Pr- 

The following summarizes all of the limits placed on the controls. 

Pi = (h — hr) > 0 (terrain limit) 

02 = (V- Vstaii) > 0 (stall velocity limit) 

Ps = ( Cimax ~ Cl(V, h)) > 0 (maximum lift coefficient limit) 

Pa = Psmax(v,h) > 0 (specific excess power) 


Using Bryson and Ho’s [1] notation, the performance index can be a function of both the 
states at the final time, z{tf ), and the final time itself, tf . In this case it will be a Mayer 
type function and more specifically it will be a linear combination of the weight of fuel 
consumed at the final time and the final time itself. 

4>{X(t f ),tf) = W f (t f ) + nt f (6.9) 

where /i is some positive constant. By increasing the value of fi the emphasis on minimizing 
time increases; when fi = 0 only fuel consumed is being minimized. Initial values of x and 
Wf are specified while a final value is specified only for x. The final constraint on range 
can be represented as 'ip = x(tf) — Xf = 0. As in Bryson and Ho [1], a new scalar function, 
$, a linear combination of the performance index and the constraints, is defined as the 
following: 

$ = p + u'i/j ( 6 . 10 ) 

The variational Hamiltonian, H , is as follows: 

H = \ x x + X w W f (6.11) 
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where A^ and A™ are the co-states. To have an optimal solution the Hamiltonian must be 
minimized over the possible set of control values. If the constraints ($) are not active then 


the following conditions apply: 


dH 


dh 


0 

0 


We shall examine the minimization of H from a geometric view based on the hodograph. 

The costate vector A (created by the costates A^ and A^) can be plotted on the same 
Wf versus x plane as the hodograph. Since the Hamiltonian is a dot product of the costate 
vector and the state rate vector, lines of constant Hamiltonian are perpendicular to the 
costate vector as shown in Figure 6.3. The direction of increasing Hamiltonian is in the 
direction of the costate vector. 


The optimality condition requires that we find a constant H line that separates the 
hodograph plane into a part that contains the achievable state-rates and a part that has 
no achievable state-rates. This line will be Tangent 5 to the envelope of the hodograph in 
Figure 6.2. 

The transversality conditions [1] include the requirement that H = — fi. For the fuel 
optimal (only) case we have fi = 0 and this means the constant H line must go through 
the origin. 

Combining the tangency requirement and the H = 0 requirement we find that all points 
along the ‘flat part 5 of the hodograph satisfy the optimality condition. This means that in 
cruise approximation, minimum fuel flight with no winds always occurs at the minimum 
drag velocity. The central part of our NLP solutions are nearly in equilibrium and display 
this feature. 


For positive fi the separating H line must be below the origin (Wf < 0 at x = 0). This 
leads to a tangency point on the envelope of the hodograph which is also on the maximum 
altitude curve as is seen in Figure 6.4. When /x is large, the solution is at the upper bound 
of altitude. As /x is decreased from that value the constant H line moves closer to the origin 
and the solution point moves towards the linear portion of the envelope curve as shown in 
Figure 6.4. 

Some insight into why certain behavior is seen in some of the numerical results was 
gained by this analysis. First of all, it was found that the airplane should fly at the 
minimum drag velocity in the absence of winds and with no emphasis on time. In the cases 
that were run with no winds and /x = 0 the airplane did fly at or very close to this velocity. 

Another important result of this analysis was finding that the envelope of the cruise 
model hodograph is linear and that all points along that line are extremal. When fi — 0, 
there is no unique answer. The cruise model says to fly at the conditions for any point 
along that linear extremal on the hodograph. These points correspond to the points along 
a constant weight curve in Figure 6.1 which means the airplane can fly at any altitude as 
long as it flies at the minimum drag velocity. 
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Chapter 7 

Summary and Conclusions 


An energy model was developed for flight in a horizontal wind field. A Non-Linear Pro- 
gramming problem was formulated to optimize the flight path with a performance index 
based on a weighted sum of fuel-consumed and final time. Constraint equations were de- 
rived for kinematics, energy, and mass. The energy constraint included energy exchange 
due to turns made at the waypoints. 

A prominent feature of the results in Section 5.2 can be interpreted as a three dimen- 
sional version of Zermelo’s problem [1]; the winds act like currents and the vehicle flies in 
a path such that the winds can help it to its destination. 


7.1 Throttle Oscillations 

The oscillation in throttle setting seen in some of the results could possibly be due to the 
way the integral is approximated in the mass constraint equation. To approximate the 
integral a simple average was used: 

2 (^ r [^/rlmaa; H“ 

In the near cruise condition the maximum flow-rates are nearly identical so that the integral 
is essentially At[Wf\ m a X * (£r + &)/ 2. Any change in the pair £ r , that produces the same 
average will produce the same fuel flow-rate. This argument also applies to the power 
available. 

As an alternative formulation, perhaps the left value of the throttle setting, £*, should 
have been used for the whole calculation since the controls are assumed to be piecewise 
constant, not piecewise linear. 
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7.2 Effects of Energy Interchange 

The effects of the energy interchange are more subtle. The total energy change between 
waypoints (E r — Ei) is the sum of several terms: specific excess power term with no winds 
from first integral in equation (3.21), energy terms due to winds from last two integrals in 
equation (3.21), energy due to banked lift in the turn from equation (3.25), and energy due 
to maneuvering drag from equation (3.31). Typical ascending and near cruise condition 
waypoints were selected from the three dimensional problem with the analytic wind model. 
The waypoints selected from the ascending portion of the flight were 120.6 nmi apart from 
each other in the horizontal direction and were 11,040 ft apart in altitude. These waypoints 
had the following energy interchanges between them: 

E r -Ei = 11,084.6 ft 
JP 8 dt = 11,023.7 ft 
SP s Pdt- JPdh = 60.96 ft 
Banked lift term = 0.509 ft 
Maneuvering drag term = -0.00561 ft 

The turning terms were not significant (about 0.51% of total energy change between two 
waypoints) in climbing portion of the flight. 

The waypoints selected from the near level portion of the flight were 328.2 nmi apart 
from each other in the horizontal direction and were 1,030 ft apart in altitude. These 
waypoints had the following energy interchanges between them: 

E r -E x = 6.57 ft 
f P s dt = -4.237 ft 
fP s Pdt-fPdh — -3.080 ft 
Banked lift term = 5.391 ft 
Maneuvering drag term = -0.0299 ft 

The term due to the banked turn is significant (about 82% of total energy change between 
two waypoints) in this near level portion of the flight. 

The change in total energy between the two waypoints is greater in the ascending 
portion of the flight due to a greater difference in altitude and velocity. There isn’t much 
difference between the altitudes of the waypoints in the near level portion of the flight so 
the energy terms due to turning make up a greater percentage of the total energy change. 

7.3 Software 

It was seen that the computer run time generally decreased as emphasis on flight time (//) 
increased. This was due to the solution point moving away from the linear portion of the 
hodograph so that there was a more definite and specific solution as was shown in Figure 
6.4. As (i increased the vehicle went towards the maximum limits on throttle (£ = 1) and 
altitude for a greater portion of the flight. With less emphasis on time the solution was on 
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the linear portion of the hodograph so any point along that line was a solution so it took 
more iterations to come to a specific one. 

The results of this paper can be used in the flight planning of the Theseus vehicle by 
having the human operator choose from family of solutions depending on a desired landing 
time. A family of solutions such as the one in Figure 5.18 can be created in the time it 
takes to run about 4 to 5 cases of varying fi factors. This turns out to be only about 1 to 
5 minutes, so the family can be found quickly and a flight plan selected while the airplane 
is flying. 

Finally, the computer code developed for this paper can be used in general for flight 
planning of any vehicle of this class flying in a horizontal wind field. The code is modular 
and requires a subroutine for the specific aircraft engine model where the input is the 
current air velocity and altitude and the output is the maximum thrust power available 
and the maximum fuel flow rate. A subroutine for the winds is also necessary. The wind 
model is given an altitude and distance from the Pole and it returns a wind speed. Some 
specific data for the aircraft is also needed such as the parabolic drag data. 
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Appendix A 

Point-Mass Model 


Our dynamic model can be simplified for the case of ‘low’ speed, ‘short’ range flight ‘near’ 
the surface of the Earth. Therefore an approximation of a flat, non- rotating Earth is made 
so the analysis begins with the hypothesis that a reference frame with its origin at a fixed 
point on the Earth’s surface is an inertial reference frame (called the Earth-Fixed Inertial 
Frame). Newton’s Second Law implies that the equation 

F = m%g\j (A.l) 

describes the motion of the mass point at r when subjected to the net applied force F at the 
current mass, m. Our first task will be to relate the second derivative term to meaningful 
engineering quantities, such as speed and flight-path angle. We will also develop some 
useful kinematics. 


A.l Relative Motions 

To carry out this analysis in an efficient manner we consider a general situation with two 
reference frames (cf Sec. 5.1 in [14], see Fig.A.l). It may be convenient to think that Tj 
is the fixed reference frame and that Tm is the moving reference frame. The translational 
motion of Tm with respect to Tj is given by the (time- varying) position vector R(t). The 
rotational motion of Tm with respect to Tj is given by the (time- varying) angular velocity 
vector uj{t). 

Following developments in ([15], pp 383-386) we write 

f=R + p (A.2) 

and find that: 

1n\i = Tt\i + i$\M + wxp (A.3) 
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and that: 


JW I 7 — \i + d)xp + ujx(u;xp) + 2 ujx^f\M + 772 I m- (A. 4) 

Note that with R = 0 (so that r(t) = p(i)) equation (A. 3) relates the time-derivative 
operation in one frame to that in a second frame. In the following we will be applying 
these principles to the flight motion problem. It is neccesary to clearly define the various 
reference frames; in particular we must describe the appropriate c 0 relating the relative 
angular motion. We have already agreed that the Earth-Fixed-Inertial system provides our 
inertial frame and we will now introduce several other frames. 

A. 2 Reference Frames 

The analysis in the section closely follows the development in [16]. A reference frame is 
arranged with origin Q at some fixed point on the surface of the Earth, and with x axis 
pointing North, y axis pointing East, and z axis pointing in the direction of gravity. This 
is the Earth- fixed frame as mentioned earlier and the symbol T e is used to denote it. 

The next reference frame - the local horizontal frame , J~h - will have its origin at the 
mass center of the flight vehicle. The position of the vehicle can be reckoned from the 
origin G of frame J- e . Specifically, the location can be given in terms of a North/South 
position, x, an East/West position, y, and an altitude h above the Earth’s surface. The 
coordinate axes are arranged so that they are parallel to the Earth-fixed axes. 

There is no transformation needed from T e to J~h since the two frames are aligned just 
as J~i and J~m appear to be aligned in Figure A.l. 

We now recall that for the wind- axes reference frame - T w - the origin is at the mass 
center of the flight vehicle, the x w axis is aligned with the velocity of the vehicle (relative to 
the surrounding air mass), the z w axis is in the symmetry plane of the vehicle, downward 
in the normal flight attitude, while the y w axis completes the right-handed orthogonal set. 
The relative orientation of the frame T w with respect to the local horizontal frame J~h is 
described by the velocity-yaw angle x (heading angle), the velocity-pitch angle 7 (flight-path 
angle), and the velocity-roll angle fi (bank angle). Composing these elementary rotations 
we find 

cos 7 cos x cos 7 sin x — sin 7 

Lh^w = sin fi sin 7 cos x ~ cos p sin x sin fi sin 7 sin x + cos p cos x sin fi cos 7 (A. 5) 

cos fi sin 7 cos x + sin p sin x cos p sin 7 sin x ~ sin fi cos x cos fi cos 7 

The relative angular velocity between Th and T w is due to changes in the angles x, 7 
and ji. The usual analysis leads to: 

^ w/h = Pw lw + Qwjw T T'w h w j (A. 6) 
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where 


Pw " 


1 0 — sin 7 


’ a " 

q w 

= 

0 cos fi sin fi cos 7 


7 

r w _ 


0 — sin fi cos fi cos 7 


. X . 


u)hje — 0, implying that u w /h = & w /e- The angular rate of T w relative to T e is the same as 
that of T w relative to Th- 


A. 3 Kinematics 

The position of the flight vehicle has been described in terms of x (positive to the North), 
y (positive to the East) and h (altitude). Clearly, the velocity will be related to the rate 
of change of these quantities. Using the point notation from Figure A.l and letting the 
Earth- fixed frame be J-j and the horizontal frame be J~m in the figure the following is true: 


y — d (QO) I 

VI — dt \ e 

d(^ + (AS) + (SO) i 

dt \ e 

d x tfa | | d y jh | , d — h kh I 
dt l e _r dt l e ' dt \ e 

= x i h + y j h - hk h (A. 8) 

where Vi is the inertial velocity and is the vector sum of the airplane’s velocity with respect 
to the surrounding air and wind velocity so that 

Vj = V + V w . (A.9) 

In this case the winds are in the horizontal direction only so 

Vw — V W x T V W y jh (A. 10) 

V wx and V wy are the horizontal components of the wind velocity. Solving equation (A.9) 
for V and substituting in equations (A. 8) and (A. 10) yields the following equation: 

V — x %h T y jh fokfa Vwx ^ h V W y jh (A. 11) 

We also have 

V = V i w = V(cos'ycosxh + cos7sinxj^ — sin7 k h ), (A.12) 
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The two equations (A. 11 and A. 12) provide equivalent expressions for the vector V. Equat- 
ing components leads us to: 

= V cos 7 cos x - V wx 


x 


y 

k 


Vu 


wy 


= V cos 7 sin % 

= V sin 7 

The expression (A. 13) provides the kinematic relation for our problem. 


(A. 13) 


A. 4 Inertial Acceleration 


We now proceed with the main task of working out the expression for the acceleration term 
in (A.l). 

From the definition of the inertial velocity in Equation (A. 9) the inertial acceleration is 
now computed as: 

#|e = ^|e+ ^le. (A. 14) 

Before proceeding with a coordinate representation of the vector terms in (A. 14), it is 
worth considering the question of which reference frame might be best for this purpose. 
While several choices are reasonable, the most common in the study of flight paths is the 
wind-axes system. One reason for this choice is the representation of the aerodynamic force. 

Based on (A. 14) the inertial acceleration will be composed of two terms. Since we opt 
for wind-axes, it is natural to compute the first of these as: 


dV 
d t 


dV 
d t 


1“ ^ w/e X V , 

and expanding the cross-product term we find 

77 " |e — Vi w T V (r w jw 

The second (wind) term in (A. 14) is found to be: 


dV m 
d t 


dV m 
d t 


V WX %h T Vyjy Jh 


(A. 15) 


(A. 16) 


where V wx and V wy are the horizontal components of the wind acceleration. Transforming 
into the wind axes components, 


(A. 17) 


When the vector representation is transformed into the wind axis system, there can be 
a third component of the wind acceleration in this system if the airplane is banking or 
pitching. 


' W x ' 


V wx 1 

w 2 

. W 7 3 . 

— ^h\-^w 

1 

't^ ° 

i 
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A. 5 Force Representation 

We commonly decompose the net force acting on the flight vehicle into aerodynamic , propul- 
sive and gravitational contributions. 

F = A + T + W (A. 18) 


A. 5.1 Aerodynamic Force 

Symmetric flight is being assumed which means no side slip, (3 = 0, and no side force, 
Q = 0. The components of the aerodynamic force in the wind-axes systems are then 
denoted by drag and lift only; that is 

A = — ( Di w + Lk w ). (A. 19) 

The D(-) and L(-) functions are specific to the vehicle whose flight performance we wish 
to model. 

A. 5. 2 Propulsive Force 

In this case the thrust direction is always assumed to be in the aircraft velocity direction. 
Specifically, we find 

f = Ti w (A.20) 


A. 5. 3 Gravitational Force 

The gravitational force is easily represented in the J~h frame as W = mg(h)kh . In this 
model acceleration due to gravity, g , is assumed to be constant at all altitudes. The 
components in the T w frame can be found from the third column of the transformation 
matrix (A. 5) which leads to 

W = mg sin j i w + sin /x cos 7 j w + cos pi cos 7 k w ) . (A. 21) 


A. 6 Assembling the Pieces 

Our objective is a dynamic model (in this case a system of ordinary differential equations), 
that will describe the time evolution of the position/velocity of the vehicle. We have chosen 
to parameterize the position in terms of (#, y, h ) and the velocity in terms of (V, 7, y). 
The differential equations describing the evolution of (rr, y, h ) have been given explicitly 
in (A. 13). The model for the evolution of (F, 7, x) 18 more intricate. 
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Our basic equation is (A.l). The description of the acceleration term has been done in 
Section (5), while the force terms were analyzed in Section (6). We have elected a wind-axes 
resolution of the vectors. To describe the process for extracting V, 7, and x, it is helpful 
to write our version of (A.l) as: 

vn (V i w + V (t w j w — Qwk w ^ = A + T + W — m (\V\ i w + IV2 j w + W3 k w ^ . (A. 22) 

The wind- axes representations of the various terms on the right-side of (A. 22) are given in 
the previous sections. Each term is computable from the current values of 

(x, y, h, V, 7, % , m) and (/x). 

The former seven variables are called state variables - these are the variables whose evolution 
is governed by the differential equations. The latter variable, along with throttle-setting are 
called control variables - these are viewed as selectable by the pilot. The pilot (or auto-pilot) 
effects the motions of the vehicle through the choice of the (time-varying) controls. 

The i w component of (A. 22) may be solved directly to obtain V. The j w component 
yields r w , while the k w component yields q w . With these known, the bottom two rows of 
the matrix equation (A. 7) can be solved for 7 and %. The seventh state- variable is the 
mass of the vehicle - its rate of decrease is the fuel mass-flow-rate which must be described 
as part of the propulsive model. 

The last two rows of (A. 7) can be manipulated to produce: 


7 = cos fi q w — sin /x r w 
cos 7 x = sin /x q w + cos /x r w 


Combining this with (A. 22) leads to our model: 


V 

7 

X 


(T — D — mg sin 7) jm — W\ 

L • 

— cos /x — g cos 7 + W3 cos /x + W2 sin ji 
m 

— sin ll — W 2 cos a + W3 sin a /(V cos 7) 

m J 


/V 


(A.23) 


The seventh state can be written as the rate of fuel consumed since the weight of fuel varies 
in the same way the weight of the entire aircraft varies. This seventh state is 

Wf = — (A. 24) 

T] 

where Wj is the weight of fuel consumed, c is the specific fuel consumption, P a is the power 
available, and rj is the propeller efficiency. 

Equations (A. 13), (A.23), and (A. 24) make up the point-mass model for symmetric 
flight over in flat non-rotating Earth in a horizontal wind field. 
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Appendix B 

Calculation of Power Required 


A procedure similar to ARPS [11] is followed to determine the PIW-VIW curve for the 
aircraft. A PIW-VIW curve is a nondimensional power versus velocity curve which is good 
for a certain aircraft configuration at any weight and altitude. Assuming a parabolic drag 
polar as in equation (3.27) and that L = W , 


D _Cd _C Do + k C L 2 
W C L C l 


(B.l) 


After substituting Cl = j — 


into equation (B.l) where q is the dynamic pressure we get: 


D _C Do k W/S 
W W/S q+ q 


(B.2) 


Using the lift coefficient for minimum drag, CL m d = JCdo/^, the velocity for minim u m 
drag is found to be 


Vmd = 


2 (W/S) 

\ PosfCoojk 


K 


emd 


a 


(B.3) 


where the e subscript in V emd is for equivalent velocity so that V e = V and & = -£-. 
After substituting q = \p Q V e 2 and doing some algebra using V emd in equation (B.2) the 
following expressions for the drag to weight ratio are obtained: 


D_ 

w 


[{L/ D)max\ 


+ £ 


u 2 + \ 


(BA) 


where u = V/Vmd • 
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The power required to overcome the drag force is P r = D • V . Now the ^ expression 
from equation (B.4) is substituted into P r and both sides are multiplied by ^fo to get the 
following expression which is good for any altitude. 


P r \f<7 




(W V e 


emd ) 


« 3 + 


u J 


(B.5) 


The equivalent velocity for minimum drag still depends on weight so a distinction is made 
between weights. Wt will be the “test” or actual weight and W s will be a standard weight 
(say the full weight of the aircraft). In all equations until now W has been the test weight. 
A new velocity parameter, independent of weight, is now introduced. 


V 


k l 2(W a /S) 

Cdo V Po 


(B.6) 


V is the air speed for minimum drag at the standard weight and at standard sea level 
density. We have V ern d = ^ is substituted into equation (B.5) and both sides are 

divided by ( WtjW s ) 3//2 and the scaling power, W s Vto get a nondimensional expression for 
all weights and altitudes. 

(Wt/W s f /2 IWD) max ] 

Ws V) = 2 

The parameters PIW and VIW are now defined nondimensionally as: 


tr + - 


(B.7) 


PIW 

VIW 


Pr\W 

(Ws V) 

(v/v) sft 

x jWt/W , 


Note that PIW only depends on VIW. The only data needed about the aircraft in order to 
draw its PIW-VIW curve, assuming a parabolic drag polar, is (L/D) max as shown in the 
following equation obtained by substituting the definitions of PIW and VIW into equation 
(B.7). 


p jtt/ _ [(L/tymax] 

2 


VIW 3 + 


1 1 

VIW. 


(B.8) 


A value for VIW can be found by knowing the following: V/V , cr, and W t /W s . PIW 
can be calculated from VIW as in equation (B.8). Finally, the power required is obtained 
nondimensionally using values already mentioned from the following equation: 


P r = PIW (W t /W s f 2 
(Ws V) ^ 
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Appendix C 
Aircraft Data 


C.l Aircraft Dimensions 


maximum gross weight = 5511 lb 
empty weight = 3825 lb 

wing area (S) = 678 ft 2 

aspect ratio ( AR ) = 28 


C.2 Flight Characteristics 

The vehicle is to operate at high altitudes (60,000 - 80,000 feet) and low dynamic pressure, 
Q, (a few pounds per square ft) so it has been designed with a low wing loading, Cl Q = 
W/S. Since the Q is low, W/S has to be low. A parabolic drag polar, Cd = Cd q + & C\ 
is used where 

e = 0.96 
C Do = 0.0153 

k = (7rARe) -1 = 0.01184 

C.3 Engine Model 

The aircraft has two engines at 80 hp each. 

Specific Fuel Consumption (c) = 0.45 lb/hr/HP 
Propeller Efficiency = 83% 

The power available is constant up to 65,000 ft and will decrease linearly to 60% at 82,000 
ft. 
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